Supplementary Materials for XX title

I. Data preparation

## Joining, by = c("Eleve", "Groupe")

Long format

One row per outcome and per assessment time

# create long format data
dat_long <- 
  dat %>%
  pivot_longer(!c(participant, experimentatrice, eleve, group, ID, age, sexe, missing, s1_bin, s1_0_1, s2_bin, s2_0_3, s3_bin, s3_0_3, s4_bin, s4_0_4), names_to = "outcome", values_to = "value") %>%
  filter(!outcome %in% c("t0_corsi_env_nb", "t1_corsi_env_nb",
                         "t0_mdc_env_nb", "t1_mdc_env_nb"))

# identify timing of assessment
dat_long$time <- ifelse(grepl("t1", dat_long$outcome), "t1", "t0")

# remove assessment time to the name of outcomes
dat_long$outcome <- stringr::str_replace_all(dat_long$outcome, "t0_", "")
dat_long$outcome <- stringr::str_replace_all(dat_long$outcome, "t1_", "")
dat_long$outcome <- stringr::str_replace_all(dat_long$outcome, " ", "_")

dat_long <- rename_outcomes(dat_long)
dat_long$value <- ifelse(dat_long$flip == "yes",
                      - dat_long$value, dat_long$value)

Semi wide format

One row per outcome (i.e., each assessment time is stored in a different column).
We also z-scaled all the outcomes

# we translate the data + center our outcomes at t0 and t1 + assess the pre/post correlation for each of our outcomes
dat_semiwide_raw <- 
  dat_long %>% 
  pivot_wider(names_from = "time", values_from = "value") %>%
  dplyr::group_by(outcome) %>%
  dplyr::mutate(z0 = as.numeric(as.character(scale(t0))),
                z1 = as.numeric(as.character(scale(t1))),
                r = as.numeric(as.character(cor.test(~t0 + t1)$estimate))) %>% 
  dplyr::ungroup()

dat_semiwide <- rename_outcomes(dat_semiwide_raw)

Composite outcomes

Sum all z-scores by cognitive function / test

II. Table1

Demographic characteristics of our sample depending on the group.

dat_table1 <- 
  dat_comp %>%
  group_by(group) %>%
  mutate(sex = if_else(sexe == "M", 0, 1)) %>%
  summarise(mean_age = mean(age),
            mean_sex = mean(sex, na.rm = TRUE))

dat$sex <- factor(dat$sexe)
dependent = "group"

table1::table1(~ age + sex | group, data=dat, overall = "Overall",
                render.continuous=c(.="Mean / Median<br>SD<br>[min, max]"))

	Control (N=22)	Experimental (N=20)	Overall (N=42)
age
Mean / Median SD [min, max]	63.1 / 62.0 3.64 [58.3, 69.5]	63.8 / 64.2 3.33 [57.8, 69.5]	63.4 / 63.3 3.47 [57.8, 69.5]
sex
F	12 (54.5%)	11 (55.0%)	23 (54.8%)
M	10 (45.5%)	9 (45.0%)	19 (45.2%)

III. Data exploration

Outcome distribution at baseline

Inhibition

dat_inhib_dist <- subset(dat_long, grepl("Inhibition", dat_long$GEN_outcome) & time != "t1")

ggplot(dat_inhib_dist, aes(x = value)) + 
  geom_histogram(aes(y=..density..), fill = "grey", alpha = 0.8) + 
  geom_density(alpha = 0.3, aes(fill = group)) + 
  facet_wrap(. ~ outcome, scales = "free") +
  theme_bw() + 
  theme(legend.position = "bottom")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Working Memory

dat_memory_dist <- subset(dat_long, grepl("memory", dat_long$GEN_outcome) & time != "t1")

ggplot(dat_memory_dist, aes(x = value)) + 
  geom_histogram(aes(y=..density..), fill = "grey", alpha = 0.8) + 
  geom_density(alpha = 0.3, aes(fill = group)) + 
  facet_wrap(. ~ outcome, scales = "free") +
  theme_bw() + 
  theme(legend.position = "bottom")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Overal EF

dat_EF_dist <- subset(dat_long, grepl("Executive", dat_long$GEN_outcome) & time != "t1")

ggplot(dat_EF_dist, aes(x = value)) + 
  geom_histogram(aes(y=..density..), fill = "grey", alpha = 0.8) + 
  geom_density(alpha = 0.3, aes(fill = group)) + 
  facet_wrap(. ~ GEN_outcome, scales = "free") +
  theme_bw() + 
  theme(legend.position = "bottom")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Other

dat_othr_dist <- subset(dat_long, !grepl("memory", dat_long$GEN_outcome) & !grepl("Inhibition", dat_long$GEN_outcome) & !grepl("Executive", dat_long$GEN_outcome) & time != "t1")

ggplot(dat_othr_dist, aes(x = value)) + 
  geom_histogram(aes(y=..density..), fill = "grey", alpha = 0.8) + 
  geom_density(alpha = 0.3, aes(fill = group)) + 
  facet_wrap(. ~ outcome, scales = "free") +
  theme_bw() + 
  theme(legend.position = "bottom")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

IV. Main analysis

Code

Table

res_main = res_main %>% mutate(across(where(is.numeric), round, digits=3))
DT::datatable(res_main, 
              rownames = FALSE,
              extensions = "Buttons",
              options = list(  # options
                scrollX = TRUE,
                dom = c('tB'), 
                buttons = c('copy', 'csv', 'excel', 'pdf', 'print'),
                scrollY = "600px", 
                pageLength = 30,
                order = list(3, 'asc'),
                columnDefs = list(
                  list(width = '70px',
                       targets = "_all"),
                  list(className = 'dt-center', 
                       targets = "_all"))))

Forest plot

tab.main <- data.frame(
  'General outcome' = res_main$outcome,
  TOST = res_main$TOST)

metaviz::viz_forest(x = res_main[, c("d", "d_se")], 
           variant = "classic",
           col = "Greys", xlab = "SMD", annotate_CI = TRUE, 
           study_table = tab.main,
           type = "study_only",
           text_size = 2.5,
           x_limit = c(-1.5, 1.5),
           x_breaks = seq(-3, 3, 1))

Box plots

preparing data

dat_plot <- 
  dat_comp %>%
  pivot_longer(
   cols = c(z0, z1),
   names_to = "time",
   values_to = "value",
  )
dat_plot$time[dat_plot$time == "z0"] <- "t0"
dat_plot$time[dat_plot$time == "z1"] <- "t1"

Inhibition

dat_inhib <- subset(dat_plot, grepl("Inhibition", dat_plot$GEN_outcome))

ggplot(dat_inhib, aes(x = time, y = value, color = group, fill = group)) +
  # geom_line() +
  geom_boxplot(alpha = 0.2) +
  facet_grid(. ~ GEN_outcome,  scales = 'free') + 
theme_bw() + 
  theme(legend.position="bottom")

Working memory

dat_memory <- subset(dat_plot, grepl("memory", dat_plot$GEN_outcome))

ggplot(dat_memory, aes(x = time, y = value, color = group, fill = group)) +
  geom_boxplot(alpha = 0.2) +
  facet_grid(. ~ GEN_outcome,  scales = 'free') + 
theme_bw()  + 
  theme(legend.position="bottom")

Overall EF

dat_EF <- subset(dat_plot, grepl("Executive", dat_plot$GEN_outcome))

ggplot(dat_EF, aes(x = time, y = value, color = group, fill = group)) +
  geom_boxplot(alpha = 0.2) +
  facet_grid(. ~ GEN_outcome,  scales = 'free') + 
theme_bw()  + 
  theme(legend.position="bottom")

Others

dat_othr <- subset(dat_plot, !grepl("memory", dat_plot$GEN_outcome) & !grepl("Inhibition", dat_plot$GEN_outcome) & !grepl("Executive", dat_plot$GEN_outcome))

ggplot(dat_othr, aes(x = time, y = value, color = group, fill = group)) +
  geom_boxplot(alpha = 0.2) +
  facet_grid(. ~ GEN_outcome,  scales = 'free') + 
  theme_bw() + 
  theme(legend.position="bottom")

V. Additional analyses

1. All outcomes

Code

row = length(unique(dat_semiwide$outcome))
res_S1 <- data.frame(outcome = rep(NA, row),
                       flip = rep(NA, row),
                       n_cont  = rep(NA, row),
                       n_exp = rep(NA, row),
                       d = rep(NA, row) ,
                       d_se = rep(NA, row) ,
                       d_low = rep(NA, row),
                       d_up = rep(NA, row),
                       d_tot = rep(NA,row),
                       B = rep(NA, row),
                       SE = rep(NA, row),
                       p = rep(NA, row))
a = 0

for (out in unique(dat_semiwide$outcome)) {
  # initialize some settings and re-initializing them as NA at each loop
  a = a+1
  mod = d = means = NA
  
  # subset dataframe ----------------
  dat_i = subset(dat_semiwide, outcome == out)
  
  ## ancova pre/post adjusted
  mod = lm(z1 ~ group + z0 + age + sexe, data = dat_i)
  means = emmeans::emmeans(mod, ~factor(group))
  # identify critical information ------------------------------------------------------------
  res_S1$outcome[a] = out
  res_S1$flip[a] = unique(dat_i$flip)
  res_S1$n_cont[a] = nrow(subset(dat_i, group == "Control"))
  res_S1$n_exp[a] = nrow(subset(dat_i, group == "Experimental"))
  res_S1$B[a] = summary(mod)$coefficients[2,1]
  res_S1$SE[a] = summary(mod)$coefficients[2,2]
  res_S1$p[a] = summary(mod)$coefficients[2,4]

  # extract results ---------
  d = smd(x = means, df = dat_i, level = 95)
  
  res_S1$d[a] = d[[1]]
  res_S1$d_low[a] = d[[2]]
  res_S1$d_up[a] = d[[3]]
  res_S1$d_tot[a] = d[[4]]
}

## Warning in qnorm(1 - cer): NaNs produced

## Warning in qnorm(1 - cer): NaNs produced

## Warning in qnorm(1 - cer): NaNs produced

## Warning in qnorm(1 - cer): NaNs produced

res_S1$d_se = (res_S1$d_up - res_S1$d_low) / (2*1.96)
### TOST analysis
res_S1$TOST <- NA
for (i in 1:nrow(res_S1)) {
  res_S1$TOST[i] <- max(
    abs(tsum_TOST(m1 = res_S1$d[i], m2 = 0, sd1 = 1, sd2 = 1, 
              n1 = res_S1$n_cont[i], n2 = res_S1$n_exp[i], 
              low_eqbound = -0.8, high_eqbound = 0.8, alpha = 0.05,
              var.equal = FALSE, bias_correction = FALSE,
              eqbound_type = "raw")$smd$dhigh),
    abs(tsum_TOST(m1 = res_S1$d[i], m2 = 0, sd1 = 1, sd2 = 1, 
              n1 = res_S1$n_cont[i], n2 = res_S1$n_exp[i], 
              low_eqbound = -0.8, high_eqbound = 0.8, alpha = 0.05,
              var.equal = FALSE, bias_correction = FALSE,
              eqbound_type = "raw")$smd$dlow))
}

Table

res_S1 <- rename_outcomes(res_S1)

res_S1 = res_S1 %>% mutate(across(where(is.numeric), round, digits=3))

DT::datatable(res_S1, 
              rownames = FALSE,
              extensions = "Buttons",
              options = list(  # options
                scrollX = TRUE,
                dom = c('tB'), 
                buttons = c('copy', 'csv', 'excel', 'pdf', 'print'),
                scrollY = "600px", 
                pageLength = 30,
                order = list(3, 'asc'),
                autoWidth = TRUE,
                columnDefs = list(
                  list(width = '70px',
                       targets = "_all"),
                  list(className = 'dt-center', 
                       targets = "_all"))))

Forest plot

tab.supp<- data.frame(
  'General outcome' = res_S1$GEN_outcome,
  Outcome = res_S1$outcome,
  TOST = res_S1$TOST)

metaviz::viz_forest(x = res_S1[, c("d", "d_se")], 
           variant = "classic",
           col = "Greys", xlab = "SMD", annotate_CI = TRUE, 
           group = res_S1$GEN_outcome,
           study_table = tab.supp,
           type = "study_only",
           text_size = 2.2,
           x_limit = c(-1.5, 1.5),
           # N = tab.prim$N,
           x_breaks = seq(-3, 3, 1))

2. No demographic covariates

Code

res_S2 <- data.frame(outcome = rep(NA, 9),
                       flip = rep(NA, 9),
                       n_cont  = rep(NA, 9),
                       n_exp = rep(NA, 9),
                       d = rep(NA, 9) ,
                       d_se = rep(NA, 9) ,
                       d_low = rep(NA, 9),
                       d_up = rep(NA, 9),
                       d_tot = rep(NA,9),
                       B = rep(NA, 9),
                       SE = rep(NA, 9),
                       p = rep(NA, 9))
a = 0

for (out in unique(dat_comp$GEN_outcome)) {
  # initialize some settings and re-initializing them as NA at each loop
  a = a+1
  mod = d = means = NA
  
  # subset dataframe ----------------
  dat_i = subset(dat_comp, GEN_outcome == out)
  
  ## ancova pre/post adjusted
  mod = lm(z1 ~ group + z0, data = dat_i)
  means = emmeans::emmeans(mod, ~factor(group))
  # identify critical information ------------------------------------------------------------
  res_S2$outcome[a] = out
  res_S2$flip[a] = unique(dat_i$flip)
  res_S2$n_cont[a] = nrow(subset(dat_i, group == "Control"))
  res_S2$n_exp[a] = nrow(subset(dat_i, group == "Experimental"))
  res_S2$B[a] = summary(mod)$coefficients[2,1]
  res_S2$SE[a] = summary(mod)$coefficients[2,2]
  res_S2$p[a] = summary(mod)$coefficients[2,4]

  # extract results ---------
  d = smd(x = means, df = dat_i, level = 95)
  
  res_S2$d[a] = d[[1]]
  res_S2$d_low[a] = d[[2]]
  res_S2$d_up[a] = d[[3]]
  res_S2$d_tot[a] = d[[4]]
}
res_S2$d_se = (res_S2$d_up - res_S2$d_low) / (2*1.96)
### TOST analysis
res_S2$TOST <- NA
for (i in 1:nrow(res_S2)) {
  res_S2$TOST[i] <- max(
    abs(tsum_TOST(m1 = res_S2$d[i], m2 = 0, sd1 = 1, sd2 = 1, 
              n1 = res_S2$n_cont[i], n2 = res_S2$n_exp[i], 
              low_eqbound = -0.8, high_eqbound = 0.8, alpha = 0.05,
              var.equal = FALSE, bias_correction = FALSE,
              eqbound_type = "raw")$smd$dhigh),
    abs(tsum_TOST(m1 = res_S2$d[i], m2 = 0, sd1 = 1, sd2 = 1, 
              n1 = res_S2$n_cont[i], n2 = res_S2$n_exp[i], 
              low_eqbound = -0.8, high_eqbound = 0.8, alpha = 0.05,
              var.equal = FALSE, bias_correction = FALSE,
              eqbound_type = "raw")$smd$dlow))
}

Table

res_S2 = res_S2 %>% mutate(across(where(is.numeric), round, digits=3))
DT::datatable(res_S2, 
              rownames = FALSE,
              extensions = "Buttons",
              options = list(  # options
                scrollX = TRUE,
                dom = c('tB'), 
                buttons = c('copy', 'csv', 'excel', 'pdf', 'print'),
                scrollY = "600px", 
                pageLength = 30,
                order = list(3, 'asc'),
                columnDefs = list(
                  list(width = '70px',
                       targets = "_all"),
                  list(className = 'dt-center', 
                       targets = "_all"))))

Forest plot

tab.S2 <- data.frame(
  'General outcome' = res_S2$outcome,
  TOST = res_S2$TOST)

metaviz::viz_forest(x = res_S2[, c("d", "d_se")], 
           variant = "classic",
           col = "Greys", xlab = "SMD", annotate_CI = TRUE, 
           study_table = tab.S2,
           type = "study_only",
           text_size = 2.5,
           x_limit = c(-1.5, 1.5),
           x_breaks = seq(-3, 3, 1))

3. Beebot mastering

dat_eval = subset(dat, s4_bin == 1)
paste0(round(nrow(dat_eval)/sum(!is.na(dat$s4_bin))*100, 3), "% of the sample succeed at a beebot mastering test after the last training session, with a mean score of ", round(mean(dat_eval$s4_0_4), 3), "/4")

## [1] "89.474% of the sample succeed at a beebot mastering test after the last training session, with a mean score of 3.824/4"

dat_semiwide$z_change = dat_semiwide$z1 - dat_semiwide$z0
dat_s3 = subset(dat_semiwide, !is.na(s4_0_4))
dat_s3$s4_bin = factor(dat_s3$s4_bin)

summary(lmerTest::lmer(z_change ~ ordered(s4_0_4) + (1|participant) + (1|outcome), data = dat_s3))

## boundary (singular) fit: see ?isSingular

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: z_change ~ ordered(s4_0_4) + (1 | participant) + (1 | outcome)
##    Data: dat_s3
## 
## REML criterion at convergence: 967.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1003 -0.4310 -0.0649  0.4627  3.6076 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  participant (Intercept) 0.01782  0.1335  
##  outcome     (Intercept) 0.00000  0.0000  
##  Residual                1.14216  1.0687  
## Number of obs: 323, groups:  participant, 19; outcome, 17
## 
## Fixed effects:
##                   Estimate Std. Error       df t value Pr(>|t|)
## (Intercept)        0.14038    0.11677 15.00000   1.202    0.248
## ordered(s4_0_4).L -0.15088    0.21720 15.00000  -0.695    0.498
## ordered(s4_0_4).Q -0.06353    0.23355 15.00000  -0.272    0.789
## ordered(s4_0_4).C  0.13878    0.24882 15.00000   0.558    0.585
## 
## Correlation of Fixed Effects:
##             (Intr) o(4_0_4).L o(4_0_4).Q
## or(4_0_4).L -0.618                      
## or(4_0_4).Q -0.169 -0.431               
## or(4_0_4).C  0.093 -0.102     -0.398    
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular

res_s3_bin = lmerTest::lmer(z_change ~ s4_bin + (1|participant) + (1|outcome), data = dat_s3)

## boundary (singular) fit: see ?isSingular

emmeans(res_s3_bin, ~s4_bin)

##  s4_bin emmean     SE    df lower.CL upper.CL
##  0      0.3562 0.1845 16.80  -0.0335    0.746
##  1      0.0254 0.0633  9.07  -0.1176    0.168
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

ggplot(dat_s3, aes(x = s4_0_4, y = z_change)) + 
  geom_jitter(width=0.2, size = 2, alpha = 0.4) + 
  theme_bw()

ggplot(dat_s3, aes(x = factor(s4_bin), y = z_change)) + 
  geom_jitter(width=0.1, size = 2, alpha = 0.4) + 
  theme_bw()

4. Individual trajectory

dat_agg <-
  dat_long %>%
  # pivot_wider(names_from = "time", values_from = "value") %>%
  dplyr::group_by(outcome) %>%
  dplyr::mutate(z = scale(value)) %>%
  dplyr::ungroup() %>%
  dplyr::group_by(time, ID, group) %>%
  summarise(value_z = sum(z))

## `summarise()` has grouped output by 'time', 'ID'. You can override using the
## `.groups` argument.

ID <- 
  dat_agg %>%
  group_by(ID) %>%
  arrange(time) %>%
  mutate(diff = value_z - lag(value_z, default = first(value_z))) %>%
  filter(time == "t1") %>%
  mutate(Improvement = if_else(diff > 0, "Improvement", "Deterioration"))%>%
  dplyr::select(ID, Improvement) %>%
  filter(!is.na(Improvement))
 
dat_line <- merge(dat_agg, ID, by = "ID")

Plot

ggplot(dat_line, aes(x = time, y = value_z, color = group, group = ID)) +
  geom_point(size = 2.5, alpha = 0.4) + 
  geom_line(size = 1, linetype  = "dotted", alpha = 0.8) +
  theme_bw() + 
  facet_grid(group ~  Improvement) +
  theme(legend.position = "none")